function dS  = CR3BP_2D(~,S,u)
    x = S(1);
    y = S(2);
    vx = S(3);
    vy = S(4);
    %[] saperate states into each components

    r1 = ((x+u)^2+y^2)^(1/2);
    r2 = ((x+u-1)^2+y^2)^(1/2);
    %[] Range of satellite wrt first and second primary bodies.
    
    n = length(S);
    %[] Length of the state vector!
    
    dS = zeros(n,1);   %[] Pre allocate memory

    dS(1:2) = S(3:4);
    %[] Derivative of the position are the velosity at that time.
    
    dS(3) = x  -  (1-u)*(x+u)/r1^3  -  u*(x+u-1)/r2^3 + 2*vy;
    dS(4) = y  -  (1-u)*y/r1^3 - u*y/r2^3 - 2*vx;
    %[] Acceleration in the rotating CRTBP frame.

end